;pro staralign1,image_stack,resolution=resolution,maxshift=maxshift,$
 ;   hmin=hmin,fwhm=fwhm, roundlim=roundlim, sharplim=sharplim, atv=atv,$
  ;  median=median, offsets=offsets, xoffhist=xoffhist, yoffhist=yoffhist,$
  ;  _extra=extra


function staralign1,image_stack,resolution=resolution,maxshift=maxshift,$
    hmin=hmin,fwhm=fwhm, roundlim=roundlim, sharplim=sharplim, atv=atv,$
    median=median, offsets=offsets, xoffhist=xoffhist, yoffhist=yoffhist,$
    _extra=extra
;+
; NAME:
;   STARALIGN1
; PURPOSE:
;       Align a stack of images, via shifts, based on locations
;       of stars in each frame.  Use on frames which are within several pixels
;       of aligned to achieve sub-pixel registration.
; EXPLANATION:
;       Locates the point sources in each frame of an
;       image stack.  For each frame (except first), matches point source
;       locations to those of first frame and computes offsets for each.
;       The X and Y star offsets that occur most frequently are taken as the
;       optimal image offsets.  Each image (except the first) is
;     interpolated on new, shifted pixel centers to remove these offsets
;
; CALLING SEQUENCE:
;       Result = FINDOFFS( image_stack [, RESOLUTION=, MAXSHIFT=, HMIN=,
;                         FWHM=, ROUNDLIM=, SHARPLIM=, OFFSETS=, XOFFHIST=,
;                         YOFFHIST=, ATV=, /MEDIAN] )
;
; OUTPUT:
;     Result is an image stack of the same dimensions as the original stack,
;     where each frame has been interplated from the corresponding original
;       frame, such that the star locations match those of the first frame.
;
; INPUTS:
;       image_stack = a 3-D array, a stack of 2-D images with IMAGE_STACK[*,*,0]
;       containing the reference frame against which the others will be shifted.
;
; OPTIONAL INPUT KEYWORDS:
;       RESOLUTION = the sub-pixel resolution desired for the registration,
;                    in pixels; default = 0.1 .  If set too fine, the
;                    star offsets histogram will be noisy and the correct offset
;                    may not be properly chosen.
;       MAXSHIFT = the maximum pixel offset between images; default = 5.0
;
;       HMIN = the threshold in counts above which stars will be sought by the
;              FIND routine; default = 18000.
;       FWHM = the characteristic width in pixels of stars in the images;
;              default = 3.0
;       ROUNDLIM = the allowed limits on "roundness" for stars, as defined in
;                  FIND.PRO; default = [-1.0,1.0]
;       SHARPLIM = the allowed limits on "sharpness" for stars, as defined in
;                  FIND.PRO; default = [0.2,1.0]
;       MEDIAN  -- set this keyword to use the median X and Y star offsets as
;                  the optimal offsets.  If not set, the location of the
;                  peak of the offsets histogram in each dimension is used as
;                  the optimal X and Y offset.
;       ATV     -- set this keyword to the number of a frame (1 <= ATV <=N,
;                  where N is the number of images in the stack)
;                  to call ATV.PRO to display the reference
;                  image (frame 0), with circles drawn on in two colors to show
;                  the locations of sources found in frame 0 and the
;                  corresponding locations from the selected frame.
;
; OPTIONAL OUTPUT KEYWORD:
;       OFFSETS = the optimal image offsets.  A 2xN array, where N is the number
;                 of images stacked in IMAGE_STACK.  OFFSETS[0,*] contains the
;                 X offsets, and OFFSETS[1,*] contains the Y offsets.
;                 OFFSETS[0,*]=[0,0] since the first image is not shifted.
;       XOFFHIST = The histogram of X star offsets for each pair of
;                 compared images--provides a useful check of the reliability
;                 of the optimal offset.
;                 A Mx2x(N-1) array, where N is the number of images stacked
;                 in IMAGE_STACK and M is the number of offset bins required to
;                 cover the range [-MAXSHIFT,MAXSHIFT] at the chosen RESOLUTION.
;                 OFFHIST[*,0,n] contains the bin center locations and
;                 OFFHIST[*,1,n] holds the histogram values for the star offsets
;                 of image (n+1) compared to the first frame (image 0).
;       YOFFHIST = The histogram of Y star offsets for each pair of
;                 compared images.
;
; PROCEDURE:
;       Call FIND to get star positions for 1st image.  Loop over all other
;       images; for each, 1) FIND star positions; 2) use SRCOR to match stars
;       between current and 1st images; 3) use list of x and y offsets for
;       matched stars to determine optimal offset (either median or peak of
;       histogram for star offsets); 4) call INTERPOLATE to shift images by
;       fractional pixels.  Locations in the new images which are outside the
;       original images are filled with NaN (Not a Number) values.
;
; NOTES:
;     additional keyword values can be passed to ATV
;       via keyword inheritance.
;
; EXAMPLE:
;     If one begins with IMAGE1 and IMAGE2, which are a pair of frames
;     aligned to within 4 pixels:
;       > image_stack = [[[image1]],[[image2]]]   ;create stack of images
;       > aligned_stack = staralign(image_stack, resolution=0.1, hmin=10000.,$
;                                   maxshift=4.0, fwhm=2.5, xoffhist=xoffhist,$
;                                   yoffhist=yoffhist, offsets=offsets)
;     ALIGNED_STACK contains the same two frames, with the second interpolated
;     so that the mode of stellar offsets is zero.  The histogram used to
;     determine the mode for this second frame can be viewed as follows:
;       > plot,xoffhist[*,0,0],xoffhist[*,1,0],psym=10
;       > oplot,yoffhist[*,0,0],yoffhist[*,1,0],psym=10,line=2
;       > print,offsets[*,1]    ;print optimal offsets for frame 1 vs. frame 0
;-
On_error,2

;   check parameters.
if n_elements(resolution) eq 0 then resolution=0.1
if n_elements(maxshift) eq 0 then maxshift=5.0
if n_elements(hmin) eq 0 then hmin=18000.
if n_elements(fwhm) eq 0 then fwhm=3.0
if n_elements(roundlim) ne 2 then roundlim=[-1.0,1.0]
if n_elements(sharplim) ne 2 then sharplim=[0.2,1.0]
if n_elements(atv) eq 0 then atv=0

; check if image stack is provided
if size(image_stack,/TNAME) eq 'UNDEFINED' then begin
    ; get a list of files to process
    image_stack=dialog_pickfile(/multi,filter='*.fits', $
       title='FITS images to align')
    if image_stack[0] eq '' then print, 'no images'
endif

; check if input is a string array
if size(image_stack,/TNAME) eq 'STRING' then begin
    imlist=image_stack
    nims=n_elements(imlist)
    ; get the array size from the first image
    img = mrdfits(image_stack[0],0,h1,/unsigned)
    imdims=size(img,/dimensions)
    ncol=imdims[0]
    nrow=imdims[1]
    tempimg=fltarr(ncol,nrow,nims,/NOZERO)
    tempimg[0,0,0]=temporary(img)
    for fn=1,nims-1 do tempimg[0,0,fn]=mrdfits(images[fn],0,h1,/unsigned)
    image_stack=temporary(tempimg)
endif

;get size of input array
ndims=size(image_stack,/n_dimensions)
if ndims ne 3 then message, 'IMAGE_STACK must be 3-D stack of 2-D images'
dims=size(image_stack,/dimensions)
nstack=dims[2]
if nstack lt 2 then $
    message, 'IMAGE_STACK must be a 3-D array of at least 2 images'

if atv ge nstack then atv=0

;call FIND to get pixel locations of stars in frame 0
find,image_stack[*,*,0],x0,y0,blah1,blah2,blah3,hmin,fwhm,roundlim,sharplim,$
    /silent
offsets=fltarr(2,nstack)
newstack=fltarr(dims)
newstack[0,0,0]=image_stack[*,*,0]   ;first frame of new stack is same as old
mincenter=float(round(-maxshift/resolution))    ;center of lowest bin
histmin=(mincenter-0.5)*resolution  ;minimum of lowest histogram bin
nbins=ceil(maxshift/resolution)*2   ;number of bins to span +/-maxshift
xoffhist=(yoffhist=fltarr(nbins,2,nstack-1))
abscissa=(findgen(nbins)-float(nbins)/2.)*resolution ;hist bin locations
for i=1,nstack-1 do begin   ;loop over all images after 0
    find,image_stack[*,*,i],x,y,blah1,blah2,blah3,hmin,fwhm,roundlim,sharplim,$
       /silent     ;find pixel locations of stars in current image
    srcor,x0,y0,x,y,maxshift,ind1,ind2,option=1       ;match stars between images
    xoffs=x0[ind1]-x[ind2] ;get offsets of all star matches
    yoffs=y0[ind1]-y[ind2]
    xh=histogram(xoffs,min=histmin,binsize=resolution,nbins=nbins) ;build hist
    yh=histogram(yoffs,min=histmin,binsize=resolution,nbins=nbins)
    xoffhist[*,*,i-1]=[[abscissa],[xh]]    ;build hist array for all images
    yoffhist[*,*,i-1]=[[abscissa],[yh]]
    if keyword_set(median) then begin
       ; use median of star offsets as optimal
       xmoff=median(xoffs)
       ymoff=median(yoffs)
       xfact=float(round(xmoff/resolution))
       yfact=float(round(ymoff/resolution))
       xfoff=xfact*resolution
       yfoff=yfact*resolution
    endif else begin
       ; take location of peak histogram value as optimal offset
       xfoff=float((where(xh eq max(xh)))+0.5)*resolution+histmin
       yfoff=float((where(yh eq max(yh)))+0.5)*resolution+histmin
       if n_elements(xfoff) gt 1 then begin ; if multiple peaks
         print,'WARNING:  '+strtrim(n_elements(xfoff),2)+$
          ' maxima in X correlation histogram:'
         print, xfoff
         xfoff=median(xfoff,/even)  ; take median of peaks as optimal offset
         print,'using median of these maxima.  medain = '+strtrim(xfoff,2)
       endif else begin
         xfoff=xfoff[0]
       endelse
       if n_elements(yfoff) gt 1 then begin
         print,'WARNING:  '+strtrim(n_elements(yfoff),2)+$
          ' maxima in Y correlation histogram:'
         print,yfoff
         yfoff=median(yfoff,/even)
         print,'using median of these maxima.  median = '+strtrim(yfoff,2)
       endif else begin
         yfoff=yfoff[0]
       endelse
    endelse
    xlocs=findgen(dims[0])-xfoff   ;locations for interpolation
    ylocs=findgen(dims[1])-yfoff
    offsets[*,i]=[xfoff,yfoff]
    ; interpolate on grid of (x,y) values
    arr=interpolate(image_stack[*,*,i], $
       xlocs,ylocs,cubic=-0.5,missing=!values.f_nan,/grid)
    newstack[0,0,i]=arr    ;add interpolated image to new image stack
    if atv eq i then begin ;display offsets for this frame with ATV?
       atv,image_stack[*,*,0],_extra=extra
       atvcircle,4.5,x0[ind1],y0[ind1],/scaled,color='green'
       atvcircle,4.5,x[ind2],y[ind2],/scaled,color='red'
    endif
endfor
return,newstack
end
